// Write resolutions
lc_normal     = 25000.000000;
lc_top        = 10000.000000;
lc_interface  = 1100.000000;
lc_litho      = 10000.000000;
lc_left_corner= 100.000000;

// Write subduction interface points
Point(1) = {0.000000,-0.000000,0.000000,lc_left_corner};
Point(2) = {12000.000000,-500.000000,0.000000,lc_left_corner};
Point(3) = {24000.000000,-1200.000000,0.000000,lc_left_corner};
Point(4) = {36000.000000,-2100.000000,0.000000,lc_interface};
Point(5) = {48000.000000,-3200.000000,0.000000,lc_interface};
Point(6) = {60000.000000,-4500.000000,0.000000,lc_interface};
Point(7) = {72000.000000,-6000.000000,0.000000,lc_interface};
Point(8) = {84000.000000,-7700.000000,0.000000,lc_interface};
Point(9) = {96000.000000,-9600.000000,0.000000,lc_interface};
Point(10) = {108000.000000,-11700.000000,0.000000,lc_interface};
Point(11) = {120000.000000,-14000.000000,0.000000,lc_interface};
Point(12) = {132000.000000,-16500.000000,0.000000,lc_interface};
Point(13) = {144000.000000,-19200.000000,0.000000,lc_interface};
Point(14) = {156000.000000,-22100.000000,0.000000,lc_interface};
Point(15) = {168000.000000,-25200.000000,0.000000,lc_interface};
Point(16) = {180000.000000,-28500.000000,0.000000,lc_interface};
Point(17) = {192000.000000,-32000.000000,0.000000,lc_interface};
Point(18) = {204000.000000,-35700.000000,0.000000,lc_interface};
Point(19) = {216000.000000,-39600.000000,0.000000,lc_interface};
Point(20) = {228000.000000,-43700.000000,0.000000,lc_interface};
Point(21) = {240000.000000,-48000.000000,0.000000,lc_interface};
Point(22) = {252000.000000,-52500.000000,0.000000,lc_interface};
Point(23) = {264000.000000,-57200.000000,0.000000,lc_interface};
Point(24) = {276000.000000,-62100.000000,0.000000,lc_interface};
Point(25) = {288000.000000,-67200.000000,0.000000,lc_interface};
Point(26) = {300000.000000,-72500.000000,0.000000,lc_interface};
Point(27) = {312000.000000,-78000.000000,0.000000,lc_interface};
Point(28) = {324000.000000,-83700.000000,0.000000,lc_interface};
Point(29) = {336000.000000,-89600.000000,0.000000,lc_interface};
Point(30) = {348000.000000,-95700.000000,0.000000,lc_interface};
Point(31) = {360000.000000,-102000.000000,0.000000,lc_interface};
Point(32) = {372000.000000,-108500.000000,0.000000,lc_interface};
Point(33) = {384000.000000,-115200.000000,0.000000,lc_interface};
Point(34) = {396000.000000,-122100.000000,0.000000,lc_interface};
Point(35) = {408000.000000,-129200.000000,0.000000,lc_interface};
Point(36) = {420000.000000,-136500.000000,0.000000,lc_interface};
Point(37) = {432000.000000,-144000.000000,0.000000,lc_interface};
Point(38) = {444000.000000,-151700.000000,0.000000,lc_interface};
Point(39) = {456000.000000,-159600.000000,0.000000,lc_interface};
Point(40) = {468000.000000,-167700.000000,0.000000,lc_interface};
Point(41) = {480000.000000,-176000.000000,0.000000,lc_interface};
Point(42) = {492000.000000,-184300.000000,0.000000,lc_interface};
Point(43) = {504000.000000,-192600.000000,0.000000,lc_interface};
Point(44) = {516000.000000,-200900.000000,0.000000,lc_interface};
Point(45) = {528000.000000,-209200.000000,0.000000,lc_interface};
Point(46) = {540000.000000,-217500.000000,0.000000,lc_interface};
Point(47) = {552000.000000,-225800.000000,0.000000,lc_interface};
Point(48) = {564000.000000,-234100.000000,0.000000,lc_interface};
Point(49) = {576000.000000,-242400.000000,0.000000,lc_interface};
Point(50) = {588000.000000,-250700.000000,0.000000,lc_interface};

// Add anchor points for domain
Point(51) = {0.000000,-250700.000000,0.000000,lc_normal};
Point(52) = {800000.000000,-0.000000,0.000000,lc_top};
Point(53) = {800000.000000,-50000.000000,0.000000,lc_litho};
Point(54) = {800000.000000,-250700.000000,0.000000,lc_normal};

// Write channel interface points
Point(55) = {20599.978859,0.000000,0.000000,lc_left_corner};
Point(56) = {24066.519011,-202.214842,0.000000,lc_left_corner};
Point(57) = {36083.045480,-1103.454242,0.000000,lc_left_corner};
Point(58) = {48099.503719,-2204.962810,0.000000,lc_interface};
Point(59) = {60115.880700,-3506.736861,0.000000,lc_interface};
Point(60) = {72132.163720,-5008.772099,0.000000,lc_interface};
Point(61) = {84148.340453,-6711.063647,0.000000,lc_interface};
Point(62) = {96164.398987,-8613.606076,0.000000,lc_interface};
Point(63) = {108180.327869,-10716.393443,0.000000,lc_interface};
Point(64) = {120196.116135,-13019.419324,0.000000,lc_interface};
Point(65) = {132211.753347,-15522.676860,0.000000,lc_interface};
Point(66) = {144227.229616,-18226.158790,0.000000,lc_interface};
Point(67) = {156242.535625,-21129.857500,0.000000,lc_interface};
Point(68) = {168257.662651,-24233.765060,0.000000,lc_interface};
Point(69) = {180272.602572,-27537.873274,0.000000,lc_interface};
Point(70) = {192287.347886,-31042.173715,0.000000,lc_interface};
Point(71) = {204301.891705,-34746.657775,0.000000,lc_interface};
Point(72) = {216316.227766,-38651.316702,0.000000,lc_interface};
Point(73) = {228330.350425,-42756.141644,0.000000,lc_interface};
Point(74) = {240344.254649,-47061.123684,0.000000,lc_interface};
Point(75) = {248181.258702,-50000.000000,0.000000,lc_interface};
Point(76) = {252357.936012,-51566.253883,0.000000,lc_interface};
Point(77) = {264371.390676,-56271.523309,0.000000,lc_interface};
Point(78) = {276384.615385,-61176.923077,0.000000,lc_interface};
Point(79) = {288397.607438,-66282.444375,0.000000,lc_interface};
Point(80) = {300410.364677,-71588.078495,0.000000,lc_interface};
Point(81) = {312422.885465,-77093.816860,0.000000,lc_interface};
Point(82) = {324435.168661,-82799.651047,0.000000,lc_interface};
Point(83) = {336447.213595,-88705.572809,0.000000,lc_interface};
Point(84) = {348459.020052,-94811.574093,0.000000,lc_interface};
Point(85) = {360470.588235,-101117.647059,0.000000,lc_interface};
Point(86) = {372481.918750,-107623.784091,0.000000,lc_interface};
Point(87) = {384493.012572,-114329.977814,0.000000,lc_interface};
Point(88) = {396503.871026,-121236.221099,0.000000,lc_interface};
Point(89) = {408514.495755,-128342.507074,0.000000,lc_interface};
Point(90) = {420524.888703,-135648.829130,0.000000,lc_interface};
Point(91) = {432535.052081,-143155.180924,0.000000,lc_interface};
Point(92) = {444544.988351,-150861.556384,0.000000,lc_interface};
Point(93) = {456554.700196,-158767.949706,0.000000,lc_interface};
Point(94) = {468564.190505,-166874.355358,0.000000,lc_interface};
Point(95) = {480568.853535,-175177.561154,0.000000,lc_interface};
Point(96) = {492568.853535,-183477.561154,0.000000,lc_interface};
Point(97) = {504568.853535,-191777.561154,0.000000,lc_interface};
Point(98) = {516568.853535,-200077.561154,0.000000,lc_interface};
Point(99) = {528568.853535,-208377.561154,0.000000,lc_interface};
Point(100) = {540568.853535,-216677.561154,0.000000,lc_interface};
Point(101) = {552568.853535,-224977.561154,0.000000,lc_interface};
Point(102) = {564568.853535,-233277.561154,0.000000,lc_interface};
Point(103) = {576568.853535,-241577.561154,0.000000,lc_interface};
Point(104) = {589757.921746,-250700.000000,0.000000,lc_interface};

// Make lines for left box domain
Line(1) = {1,55};
Line(2) = {55,52};
Line(3) = {52,53};
Line(4) = {53,54};
Line(5) = {54,104};
Line(6) = {104,50};
Line(7) = {50,51};
Line(8) = {51,1};

// Make lines for interface
Line(9) = {1,2};
Line(10) = {2,3};
Line(11) = {3,4};
Line(12) = {4,5};
Line(13) = {5,6};
Line(14) = {6,7};
Line(15) = {7,8};
Line(16) = {8,9};
Line(17) = {9,10};
Line(18) = {10,11};
Line(19) = {11,12};
Line(20) = {12,13};
Line(21) = {13,14};
Line(22) = {14,15};
Line(23) = {15,16};
Line(24) = {16,17};
Line(25) = {17,18};
Line(26) = {18,19};
Line(27) = {19,20};
Line(28) = {20,21};
Line(29) = {21,22};
Line(30) = {22,23};
Line(31) = {23,24};
Line(32) = {24,25};
Line(33) = {25,26};
Line(34) = {26,27};
Line(35) = {27,28};
Line(36) = {28,29};
Line(37) = {29,30};
Line(38) = {30,31};
Line(39) = {31,32};
Line(40) = {32,33};
Line(41) = {33,34};
Line(42) = {34,35};
Line(43) = {35,36};
Line(44) = {36,37};
Line(45) = {37,38};
Line(46) = {38,39};
Line(47) = {39,40};
Line(48) = {40,41};
Line(49) = {41,42};
Line(50) = {42,43};
Line(51) = {43,44};
Line(52) = {44,45};
Line(53) = {45,46};
Line(54) = {46,47};
Line(55) = {47,48};
Line(56) = {48,49};
Line(57) = {49,50};

// Make lines for channel interface
Line(59) = {55,56};
Line(60) = {56,57};
Line(61) = {57,58};
Line(62) = {58,59};
Line(63) = {59,60};
Line(64) = {60,61};
Line(65) = {61,62};
Line(66) = {62,63};
Line(67) = {63,64};
Line(68) = {64,65};
Line(69) = {65,66};
Line(70) = {66,67};
Line(71) = {67,68};
Line(72) = {68,69};
Line(73) = {69,70};
Line(74) = {70,71};
Line(75) = {71,72};
Line(76) = {72,73};
Line(77) = {73,74};
Line(78) = {74,75};
Line(79) = {75,76};
Line(80) = {76,77};
Line(81) = {77,78};
Line(82) = {78,79};
Line(83) = {79,80};
Line(84) = {80,81};
Line(85) = {81,82};
Line(86) = {82,83};
Line(87) = {83,84};
Line(88) = {84,85};
Line(89) = {85,86};
Line(90) = {86,87};
Line(91) = {87,88};
Line(92) = {88,89};
Line(93) = {89,90};
Line(94) = {90,91};
Line(95) = {91,92};
Line(96) = {92,93};
Line(97) = {93,94};
Line(98) = {94,95};
Line(99) = {95,96};
Line(100) = {96,97};
Line(101) = {97,98};
Line(102) = {98,99};
Line(103) = {99,100};
Line(104) = {100,101};
Line(105) = {101,102};
Line(106) = {102,103};
Line(107) = {103,104};
Line(108) = {75,53};

// Add flag for prescribed velocities
Physical Line(101) = {9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57};

// Add flags for subduction channel
Physical Line(102) = {59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78};

Physical Line(103) = {79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107};

// Add flag for bottom overriding plate
Physical Line(104) = {108};

// Add flag for top boundary
Physical Line(111) = {1,2};

// Add flag for right boundary
Physical Line(112) = {3,4};

// Add flag for bottom boundary
Physical Line(113) = {5,6,7};

// Add flag for left boundary
Physical Line(114) = {8};

Line Loop(10) = {-8,-7,-57,-56,-55,-54,-53,-52,-51,-50,-49,-48,-47,-46,-45,-44,-43,-42,-41,-40,-39,-38,-37,-36,-35,-34,-33,-32,-31,-30,-29,-28,-27,-26,-25,-24,-23,-22,-21,-20,-19,-18,-17,-16,-15,-14,-13,-12,-11,-10,-9};

Line Loop(20) = {9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,-6,-107,-106,-105,-104,-103,-102,-101,-100,-99,-98,-97,-96,-95,-94,-93,-92,-91,-90,-89,-88,-87,-86,-85,-84,-83,-82,-81,-80,-79,-78,-77,-76,-75,-74,-73,-72,-71,-70,-69,-68,-67,-66,-65,-64,-63,-62,-61,-60,-59,-1};

Line Loop(30) = {59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,108,-3,-2};

Line Loop(40) = {79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,-5,-4,-108};

// Make plane surface
Plane Surface(1000000) = {10};
Plane Surface(2000000) = {20};
Plane Surface(3000000) = {30};
Plane Surface(4000000) = {40};

// Make physical surface
Physical Surface(0) = {1000000,2000000,3000000,4000000};

MeshAlgorithm Surface{2000000} = 3; // only initial mesh to not insert any nodes!
